In [1]:
# Para evitar los warnings de python
import warnings
warnings.filterwarnings('ignore')
%reload_ext sage
In [2]:
# Definimos el cuerpo donde vamos a trabajar
m = 4
K = GF(2)
F.<a> = GF(2^m)
In [3]:
# Creamos el anillo de polinomios
PR = PolynomialRing(F,'X')
X = PR.gen()
N = 15
L = [a^i for i in range(N)]
print(L)
g = X^3 + X + 1
y = X^2+1
In [4]:
# Imprimimos el polinomio que hemos creado
print ("Polinomio de Goppa \n\n\t g(x)= " + str(g))
In [5]:
# Declaración de variables
rango = 3
Comenzamos creando la matriz X
In [6]:
T = matrix(F,rango,rango)
for i in range(rango):
count = rango - i
for j in range(rango):
if i > j:
T[i,j]=g.list()[count]
count = count + 1
if i < j:
T[i,j] = 0
if i == j:
T[i,j] = 1
In [7]:
print ("Matriz T: ")
show(T)
Out[7]:
Definimos ahora nuestras matriz Y y Z
In [8]:
Y = matrix(F,rango,N)
for i in range(rango):
for j in range(N):
Y[i,j]=L[j]^i
show(Y)
Out[8]:
In [9]:
Z = diagonal_matrix([1/g(L[i]) for i in range(N)]) # La matriz Z sabemos es diagonal y de tamaño n x n.
show(Z)
Out[9]:
Calculamos ahora la matriz H generadora del código
In [10]:
H = T*Y*Z
show(H)
Out[10]:
Algoritmo de Patterson:
El algoritmo de Patterson es un algoritmo clásico de tiempo polinomial que permite saber si un código es decodificable.
In [11]:
# -------------------
# Funcion auxiliar que permite descomponer un polinomio en irreducibles
# -------------------
def descomponer_polinomio(p):
# El siguiente metodo permite descomponer un polinomio p en factores irreducibles p(z) = p0 (z) + z p1 (z)
# Entrada: Polinomio p
Phi1 = p.parent()
p0 = Phi1([sqrt(c) for c in p.list()[0::2]])
p1 = Phi1([sqrt(c) for c in p.list()[1::2]])
return (p0,p1)
# -------------------
# Algoritmo de euclides extendido: Obtener MCD y los s,t que lo generan.
# -------------------
def algoritmo_euclides_extendido(self, other):
delta = self.degree() #grado de polinomio 1
if other.is_zero(): # si el polinomio introducido es
ring = self.parent() #comprobamos el cuerpo en el que trabajamos
return self, R.one(), R.zero() #mcd = mismo polinomio y devuelve un uno (s) y un cero (t) en el cuerpo que trabajamos.
# mcd (a,b) = as+bt
ring = self.parent() #comprobamos el cuerpo en el que trabajamos
a = self # guardamos una copia del primer polinomio 1 (self)
b = other # guardamos una copia del segundo polinomio (other)
s = ring.one() # guardamos en s el uno del anillo
t = ring.zero() # guardamos en t el cero del anillo
resto0 = a
resto1 = b
while true:
cociente,resto_auxiliar = resto0.quo_rem(resto1) # La funcion quo_rem de Sage devuelve el cociente y el resto. Que guardamos en Q y ring.
resto0 = resto1
resto1 = resto_auxiliar
s = t
t = s - t*cociente
if resto1.degree() <= floor((delta-1)/2) and resto0.degree() <= floor((delta)/2):
break
V = (resto0-a*s)//b
coeficiente_lider = resto0.leading_coefficient() # guardamos el coeficiente lider del resto 0
return resto0/coeficiente_lider, s/coeficiente_lider, V/coeficiente_lider
# -------------------
# Funcion que calcula la inversa de un polinomio utilizando el algoritmo de euclides de Sage
# -------------------
def inversa_g(p,g):
(d,u,v) = xgcd(p,g)
return u.mod(g)
In [12]:
# -------------------
# Funcion de decodificacion de Patterson
# -------------------
def decodePatterson(y):
# Calculamos primero el vector alpha con los elementos primitivos.
alpha = vector(H*y)
# Consideramos nuestras matrices T,Y,Z definidas así como nuestro polinomio irreducible g
polinomioS = PR(0) # Inicializa como el polinomio 0 del anillo
for i in range(len(alpha)):
polinomioS = polinomioS + alpha[i]*(X^(len(alpha)-i-1)) # Lo vamos rellenando con los alpha
vector_g = descomponer_polinomio(g) # Guardamos en vector_g el par de polinomios irreducibles
w = ((vector_g[0])*inversa_g(vector_g[1],g)).mod(g)
vector_t = descomponer_polinomio(inversa_g(polinomioS,g) + X)
R = (vector_t[0]+(w)*(vector_t[1])).mod(g)
(a11,b11,c11) = algoritmo_euclides_extendido(g,R)
# Definimos el polinomio sigma
sigma = a11^2+X*(c11^2)
# Vamos comprobando uno a uno los coeficientes de sigma
# para asi determinar el conjunto de posiciones de error E - {i tal que e_i es distinto de 0}
for i in range(N):
if (sigma(a^i)==0): # error encontrado
print ("Error encontrado en la posición: " + str(i))
y[i] = y[i] + 1
return y
In [13]:
# Definimos nuestra matriz de control de paridad H de tamaño correcto (nula en estos momentos)
H_Goppa_K = matrix(K, m*H.nrows(), H.ncols())
# Y la rellenamos correctamente
for i in range (H.nrows()):
for j in range(H.ncols()):
be = bin(eval(H[i,j]._int_repr()))[2:]
be = '0'*(m-len(be))+be; be = list(be)
H_Goppa_K[m*i:m*(i+1),j]=vector(map(int,be))
print ("Matriz de control de paridad H para nuestro código") ; show(H_Goppa_K)
# Tenemos la matriz H del código Goppa, vamos a calcular la matriz G que será de dimensión (k x n)
# Para calcularla usamos GH^{T} = 0
inversa_H_Goppa_K = H_Goppa_K.right_kernel();
# Y usamos la funcion basis_matrix() que nos permite quedarnos con la matriz G.
G_Goppa = inversa_H_Goppa_K.basis_matrix();
print ("Matriz generadora G para nuestro código") ; show(G_Goppa)
Out[13]:
Out[13]:
In [14]:
## Vector sin codificar aleatorio
u = vector(K,[randint(0,1) for n in range(G_Goppa.nrows())])
print ("Vector para codificar: ") ; show(u)
## Vector codificado
c = u*G_Goppa
print ("Vector codificado: ") ; show(c)
# Añadimos un vector c con errores en las posiciones 2 y 3
e = vector(K,N)
e[2] = 1
e[3] = 1
print ("Vector e de errores:") ; show(e)
# Lo añadimos al vector que habíamos considerado
y = c + e
print ("Vector c con errores e") ; show(y)
Out[14]:
Out[14]:
Out[14]:
Out[14]:
In [15]:
decodePatterson(y)
Out[15]:
In [0]: